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Abstract 

A numerical method is presented which aUows to com- 
pute the spectrum of the Schroedinger operator for a 
particle constrained on a two dimensional flat torus un- 
der the combined action of a transverse magnetic field 
and any conservative force. The method employs a fast 
Fourier transform to accurately represent the momen- 
tum variables and takes into account the twisted bound- 
ary conditions required by the presence of the magnetic 
field. An accuracy of twelve digits is attained even with 
coarse grids. Landau levels are reproduced in the case 
of a uniform field satisfying Dirac's condition. A new 
fine structure of levels within the single Landau level is 
formed when the field has a sinusoidal component with 
period commensurable to the integer magnetic charge. 
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1. Introduction 

The quantum mechanics of a charged particle hving on a two-dimen- 
sional torus in presence of a uniform magnetic field, orthogonal to the 
surface, has been solved years ago [H [21 [3] . The degeneration of the 
ground state coincides with the flux of the magnetic field, in units 
of the elementary flux hc/e (in this paper we shall adopt units such 
that h = e/c = 1). This is a simple example of the more general 
theorem about cohomology groups for hermitian line bundles [1], known 
in the physical literature as Dirac's quantization condition: quantum 
mechanics requires that the flux of the magnetic field across a closed 
surface must be quantized. This is also known as the Weil-Souriau- 
Kostant quantization condition. 

In this paper we present a numerical algorithm which is accurate 
enough in representing the momentum variables and it respects the 
constraints posed by differential geometry. The algorithm computes 
the spectrum of the quantum particle on the torus in presence of both 
a transverse magnetic field and a scalar potential. If the potential 
vanishes and the magnetic field is uniform the algorithm reproduces 
the known spectrum, in terms of eigenvalues and degeneration, to a 
typical accuracy of twelve digits. The effect of the potential energy 
is to split the Landau Levels; this fact is at the basis of Klauder's 
formulation of path integrals in phase space [5] and our algorithm could 
be used to explore this approach to quantization theory. The case of 
a non-uniform magnetic field and the corresponding splitting pattern 
of Landau levels can be studied using our algorithm. We consider 
the case of a sinusoidal contribution to the magnetic field in the last 
section. A peculiar fine structure emerges, which is made visible by the 
accuracy of the algorithm. This fine-structure within each Landau level 
could be dubbed Landau-Mathieu levels and it manifests itself when 
the number of oscillations of the perturbed field is commensurate to 
the quantized magnetic flux. This fact suggests that an undulating 
stationary magnetic field could be used to tune the number of states 
in the fine structure of Landau-Mathieu sub levels. 



2. The model 

Quantum mechanics on a compact surface, in the presence of a mag- 
netic field transverse to the surface, requires the introduction of either 
a singular magnetic potential (Dirac's string) or a collection of local 
potentials Aa, one for each local chart of a given atlas on the surface. 
The description in terms of local potentials is preferable for its math- 
ematical rigor [6]. The implementation of the local description within 
a numerical approach should be easily achieved in terms of finite ele- 
ments methods. In this paper we take an alternative route, working on 
a single chart, but imposing the correct (twisted) boundary conditions 
to the wave function, as we explain in the next section. 
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3. Local charts and twisted boundary conditions 

Let the torus be identified with the two-dimensional plane modulo 
the discrete subgroup of translations generated by x ^ x + Li,y ^ 
y + L2. We cover the torus with four charts defined as follows 



(1) 



^ \0<x < Li ^ \6i<x<Li + Si 

^ (0<x <Li ^ (Si<x<Li + 5i 

^ ' ^ 62 <y < L2 + 52 ^ ' \ S2 <y < L2 + S2 



In each chart we define a local magnetic potential by 
(2) yt:A, = {-lBy,^^Bx) 

(remember we use units where e/c = 1). All local potentials are defined 
in the same way, but their values are different. Within the overlaps 
of the local charts wc easily find the transition functions realizing the 
gauge transformations from one description to another. For instance, 
the chart (3 overlaps a in two distinct regions, = {61 < x^ = xp < 
Li} and X^^l = < x^ < 5i,Li < xp < Li + 5i}. In the overlap X^^^ 
the value of the potentials coincide, while in X)^p we have 

.3. Ap^Aa+ (0, \BL,) 

with Xap = \BLiy. The other transition functions are determined 

similarly. For instance in Xa-^ — {0 < ya < S2, L2 < y-y < L2 + 82} it 
holds 

A = '^« + (-1^^2,0) 

with Xaj = —\BL2X. 

Now, to build the Hamiltonian operator, which is formally given by 
the usual minimal coupling, one has to establish the transition func- 
tions proper to the local wave functions. As it is well-known these are 
obtained by exponentiating the transition functions, i.e. 

(5) i^p{x, y) = e'^«'^Va(a;, y) on 

Now take a sequence of points Si converging to {Li,y) from the left and 
a second sequence S2 converging from the right to the same point. On Si 
we have ipa = i'p ^ i'ai.Li^y)] on S2 we have — >• 4>aiS^-i v) exp{ii BLi y}. 
By continuity of 'ipp wc get a condition on V'a namely 

(6) Va(^i,l/)=e^*^^^^V^«(0,i/). 
By a similar argument we find a second condition 

(7) V'«(a;,L2) =e- ^'^^^^^^x, 0) . 

At this point we are allowed to work on a single local chart (let's choose 
Co) and the Hamiltonian is defined by 

(8) H = U-id, + \Byf + U-idy - \Bxf + V{x,y) 
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on a domain of differentiable functions satisfying Eq.s (|6]l71) as boundary 
conditions. Notice that the b.c. are only consistent if Dirac's condition 
is satisfied. To see this, compute ipi^Li, L2) by applying the b.c. in two 
different orders: 

(9) ^{Li, L2) = V(0, L2) e ¥bl.L2 ^ ^^0, 0) e ¥bl.L2 

(10) ip{Li,L2) = 7/^(Li,0)e- ^^-^^^^2 = V^(0,0)e- 

hence exp{i B Li L2} = 1. All this is well-known, but it was recalled 
here to introduce the main idea behind the algorithm we describe in 
the next section. 

4. The algorithm 

A simple code, based on a discrete approximation of partial deriva- 
tives, is easily produced; the twisted boundary conditions Eq.s (0 ED 
are implemented without difficulty. However this methods has seri- 
ous limitations in attaining good accuracies. A test run with B = 27r, 
Li = L2 = 2 performed with a 64 x 64 grid in configuration space yields 
the low energy spectrum (first 20 eigenvalues) with an average error of 
1.5%. In particular the first four eigenvalues, which should coincide 
with 71, turn out to be vr x (0.9997, 1.0082, 1.0082, 1.0419). With a finer 
mesh (128 x 128) the error improves (0.5%) but the computing time 
grows considerably (from 25 sec to ~ 400 sec). This fact encourages 
to design an algorithm with a better accuracy on partial derivatives. 
This is achieved by using a "spectral method" based on the Fourier 
transform. 

4.1. The spectral method. A very accurate representation of partial 
derivatives can be obtained by using Fourier transform, in one of its 
efficient implementations as a numerical code; we shall use FFTW [S], 
which is now included in Matlab. However, Fourier transform assumes 
a periodic wave-function, which is not the case with our problem. The 
way out is to apply FFT separately along x and y; the x transform is 
applied to the function = exp{ — ^iBxy}ip, which turns out to be 
periodic in x with period Li. The minimal coupling is then recovered 
by realizing that 

(11) i-td, + ^By)ij = ey''^y{-td,<j)) + Bytlj . 

Now the partial derivative can be computed in x— Fourier space. Simi- 
larly (f) = exp{^iB X y} is periodic in y with the right period, and we 
may compute 

(12) {-idy -\Bx)^^ e-^'^'^y^-idycj)) -Bxtfj . 

The idea is used to compute with high accuracy the action of the Hamil- 
tonian on any function satisfying the twisted b.c; this is then used as 
the unique piece of information needed by the Arnoldi algorithm to get 
the spectrum. We also have to choose an initial vector, if do not feel 
easy about a random initial vector. A function satisfying the boundary 
conditions can be constructed as follows. Choose any ipQ{x,y), e.g. a 
Gaussian centered in the middle of the rectangle of sides Li,L2. Let 
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BL1L2 = 2tcN; then the following equation defines a "good" wave- 
function: 

(13) ^(x,y)= 5^(-l)'^^"^^e^^^^^^-^^^^^^^o(a; + niLi,y + n2L2) 

ni,n2 

The series can be truncated if ipo is a Gaussian with a width small with 
respect to Lj. 

These are the ingredients which can be used to make a call to Mat- 
lab's routine eigqj, which provides a very friendly interface to the 
Arnoldi package Arpack[9]. The result is rather spectacular as we re- 
port next. 

4.2. Test runs and error estimates. We apply the algorithm to a 
grid n X n, starting with very coarse grids. In Tab.l we report the 
average error and the timings to compute the first 20 eigenvalues with 
the same data as before. 



n 


Relative Error 


Time (sec) 


8 


1.8 X 10"^ 


0.15 


10 


2.5 X 10-^ 


0.25 


12 


5.3 X 10-'' 


0.35 


16 


1.0 X 10-^'^ 


0.60 


24 


1.8 X 10-1^ 


1.35 


32 


2.4 X 10-^^ 


2.85 


64 


6.0 X 10"^^ 


24.7 



Table 1. 



As we see, the algorithm reproduces the correct spectrum (includ- 
ing degeneracy) already at very low n. The relative error saturates 
around 10~^^ which seems to be inherent to the Arnoldi algorithm as 
implemented in Matlab (routine eigs). 

In Fig.l we see a typical spectrum obtained with the algorithm. The 
degeneracy of the eigenvalues is within 10~^^, obtained with a 32 x 32 
grid. 

Let us notice that if we plug a value of B which does not respect 
Dirac'c condition, the degeneracy is broken; this fact can be interpreted 
as due to the fact that there is a spurious singular contribution to the 
magnetic field at the boundary of the local chart which breaks the 
original symmetry. 

Another check for accuracy can be performed by adding a potential 
energy ^uj'^{x'^ + y"^), in which case the spectrum is known in the limit 
of large Li and L2. In the case i? = 2,0; = 1, we get the spectrum 
E = (ni -|- ^)uji + (n2 + ^)uj2 with a relative error of 10~^^ on a 64 x 64 
grid. The presence of a potential energy requires a relatively finer mesh. 



^The Matlab code can be found at the author's web site 
http : //www ■ f is . unipr . it/ ~enrico . onof ri[ 
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Figure 1. The Landau levels with = ^BLiLo = 8, 
in units of the Larmor frequency. 



5. Fine structure Landau-Mathieu levels 

Having an algorithm which allows for accurate eigenvalue computa- 
tions is like having a microscope with higher resolution power: you can 
resolve details which would otherwise be invisible. It came then as a 
surprise, using the new algorithm, to discover a structure in Landau 
levels when the uniform magnetic field is perturbed by an undulating 
additive contribution B — > B{1 + A sin(27rz/a;/Li). Notice that bound- 
ary conditions adapted to this choice of gauge fields must be reformu- 
lated, along the hues of Sec. [31 Fig. 2 shows the splitting of the first 
Landau level which occurs at u = 4. The pattern is reproduced for 
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Figure 2. The fine structure pattern of the first Landau 
level, V = A, Nq = 16. The picture below is a blow-up of the 
portion of the rectangle in the picture above. 



other choices of parameters and it looks numerically very stable and 
degeneracy within the fine structure levels is observed numerically at 
12 digits precision. (see Fig. [3]). There are A^o = ^BLiL2 states in the 
first level; these are subdivided in finer sub-levels if Nq is a multiple 
of degeneracy is given by the greatest common divisor gcd(A'"o, 2i/), 
hence it is destroyed if Nq and u are relatively prime, but it is left 
unchanged if z/ is a multiple of A^o- 

We also explored the stability of this phenomenon with respect to 
deformation of the magnetic field, by keeping its periodicity on the 
torus, e.g. by adding a higher harmonic contribution; the pattern of 
degeneracy stays the same, only the eigenvalues are shifted (see Fig. 
4). 




Figure 3. The fine structure pattern of the first Landau 
level for = 4, TVq = 32, A = 1/10. 

The finite structure energy gap is not uniform, but a regular pattern 
emerges looking at sufficiently large Nq/u. The evidence is that the 
gaps are approximately reproduced by 

(14) En+i - En oc sin(n7r/A^o) , n + u/2 = mod {2u) , 

at least when the degeneracy pattern {u, 2z/, 2i/, 2z/, u} is realised. At 
this level, however, the study is still preliminary. 

6. Concluding remarks 

We presented a spectral algorithm which can compute the energy 
spectrum for a scalar particle on the 2-D fiat torus, subject to a trans- 
verse magnetic field and any potential energy. To realize the algorithm, 
it is crucial to implement the correct boundary conditions in order to 
be able to apply the spectral method based on Fourier transform. The 
spectrum is typically obtained to a relative error of 10~^^ even on rather 
coarse meshes. When the field deviates from uniformity in a sinusoidal 
way, we find a fine structure in the splitting of Landau levels with a reg- 
ular degeneracy pattern. The problem we considered here originated 
from the formulation of the Hamiltonian path integral introduced long 
ago by J.R. Klauder [5J; see also [7\ 
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Figure 4. The gaps in the fine structure of the first Landau 
level under a periodic deformation of the magnetic field, v = 
4, A'o = 12, a is the coupling of the higher harmonic. 
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